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Abstract 

We introduce a Predictor-Corrector type method suitable for performing many-particle Brownian Dy- 
namics simulations. Since the method goes over to the Gear's method for Molecular Dynamics simulation 
in the limit of vanishing friction, we refer to it as a Gear-like algorithm. The algorithm has been tested on a 
one-dimensional, stochastically damped harmonic oscillator model, showing that it can cover a wide range 
of friction coefficients with a high-order accuracy, excellent stability, and a very small energy drift on the 
long time scales. 

PACS numbers: 02.70.-c, 47.57.-s, 52.40.Hf, 52.25.Vy 
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I. INTRODUCTION 



Brownian Dynamics (BD) simulation method 



4|] has been widely used in studying 



problems in many dispersed systems, such as polymer solutions [|5D, colloidal suspensions BaLZLlSD, 
and, more recently, complex plasmas Ballfl,lilLll2Lll3|], and has been particularly regarded as the 



"mainstay" of colloid modeling over the past two decades. 

BD method may be regarded as a generalization of the usual Molecular Dynamics (MD) 
method, which is a viewpoint that will be held throughout the present paper. As the MD method is 
based on Newton's equations of motion, the BD method is based on their stochastic generalization, 
namely, the Langevin equation and its integral (or Kramers equation written in Hamilton form), 
d 



dt 
d 

dt 



1 



7 v + — F + A(t). 



m 



(1) 



Here, as usual, m, v and r are respectively the mass, velocity and the position of a Brownian 
particle, and F is the systematic (deterministic) force resulting from external sources and/or from 
the inter-particle interactions in a many-particle system of Brownian particles interacting with the 
ambient gas of light particles. What is different from Newton's equations is the appearance of 
dynamic friction, —7V, and the stochastic (Brownian) acceleration, A(t). As already pointed 
out by Langevin, those two terms represent complementary effects of the same microscopic phe- 
nomenon: numerous, frequent collisions between the Brownian particle and molecules in the sur- 
rounding medium. While friction represents an average effect of those collisions, the stochastic 
acceleration represents fluctuations due to the discreteness of collisions with the ambient particles, 
and is represented by a delta-correlated Gaussian white noise. When the system is in thermal 
equilibrium, the friction and the stochastic acceleration are related to the ambient temperature via 
fluctuation-dissipation theorem. Note that the friction coefficient 7 is usually regarded as constant 

w. 

The Langevin equation (OQ) may be inte grat ed numerically in a manner similar to that used for 



Newton's equations in MD simulations [|4j, 



14 



15 



1711 . and the resulting technique is usually 



called Brownian Dynamics. However, when compared to the well-established MD techniques, 
algorithms for conducting BD simulations are considerably less well developed because of their 
stochastic nature. It should be noted that, besides the term Brownian Dynamics, one finds in 
literature other names that are used to designate various techniques of numerical integration of 
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Langevin equation, such as Langevin Dynamics, or Langevin Molecular Dynamics, depending on 
the background, area of application, or simply preferences of various authors ft 1 311 . It is beyond 
the scope of this paper to discuss all such techniques and subtle differences among them in the 
derivation of the simulation formulae and/or in their implementations. Therefore, to be specific, 
we are going to stick to the definitions used by Allen and Tildesley 04|], and to the strategies 
proposed by Ermak et al. et al yj], Allen yfl, van Gunsteren and Berendsen |3J], which are so far 
the most popular BD methods used in practice. 

A straightforward method of conducting the BD simulation based on Eq. (OQ) was developed 
by Ermak et al. UJ], in which equations of motion are integrated over a time interval t under the 
assumption that the systematic force F remains approximately constant over that interval. Such 
method is often referred to as an Euler-like method in the sense that it goes over to the Euler's 
scheme of standard MD simulation, in the limit of an infinitesimally small 7. In spite of the 
low efficiency of Euler's method in the MD simulation, the Euler-like method in BD simulation 
seems to be very robust for very large 7, and is still widely used in practice because of its simplicity 
PLlZLISO- However, very small time steps are required in this method to achieve bearable inaccuracy 
and energy drift, as will be demonstrated below. 

Ermak' s method can be extended to higher order schemes by keeping higher-order terms in a 
Taylor series of the deterministic force F or, equivalently in the acceleration a = F/m, as follows 



a(t) = a(0) + a(0)t + ^a(0)t 2 



ia(0)t 3 + --- + ^a^(0)t n + -- 



(2) 



where a( n ) represents the nth-order time derivative of a. For example, when integrating Eq. (OQ), 
Allen obtained a scheme by effectively keeping up to the second-order derivative in Eq. 

0. His method is reduced to the Beeman's method for MD simulation when 7^0, and is 
therefore appropriately referred to as a Beeman-like method for BD simulation. We note that this 
method has a 3rd order accuracy for both the velocity and position, using the definition of Ref. 
111411 in which the order of an algorithm is given by the highest order of the time step in simulation 
formulae. An another algorithm for integrating Eq. £Q), which also gained popularity in practice, 
was developed by van Gunsteren and Berendsen y|], by using similar assumption as Allen yLlfi]. 
Since their algorithm is reduced to the position Verlet algorithm for MD simulation in the limit 
7 — > 0, we shall refer to it as Verlet-like method for BD simulation. It was proven P, m[ that the 
Verlet-like method is numerically equivalent to the Beeman-like method in position, whereas the 
latter has higher accuracy in velocity. 



Over the years, much effort has been invested in 



rithms for BD simulation lll8 



lialm 



21 



22, 



23 



24, 



searching for better and 



25 



26, 



27 



28 



30, 



31, 



ligher-order algo- 



3211 . For example, 



similar idea of truncating the Taylor expansion was adopted by Helfand Ill8h . Iniesta and Garcia 
de la Torre Q22Q, and Honeycutt 1124b to obtain second-order Runge-Kutta-like methods, and more 
recently by Her shkovitz 1125 [to obtained a fourth-order Runge-Kutta-like method. However, most 



of those algorithms II18 . 



22j,|24, 



25 



I, 3] are based on direct extensions of the deterministic Runge- 
Kutta algorithms to include some stochastic terms [24], or on solving the deterministic part of 
the Langevin equation using standard Runge-Kutta package |25|]. Consequently, such approaches 



require more than one evaluation of the deterministic force per time step, which clearly reduces 



14, 



16. 



12b. 



their efficiency for simulating problems in many-particle systems [|4J, 

More recently, Brarika and Heyes [12611 developed a series of algorithms to improve the effi- 
ciency of BD simulation based on a finite step-size expansion for stochastic differential equations 
(SDEs), which may be therefore also classified as second-order Runge-Kutta-like methods. It was 
found that such algorithms could significantly increase the efficiency of simulation |26J]. However, 
they are mainly designed for simulating over-damped systems, such as colloidal suspensions, as 
they are based on the so-called position Langevin equation 12611 . m which the variation of the ve- 
locity [i.e., the left-hand-side of the second line in Eq. ©] is neglected IU8I1 . Such simplification 
appears to be quite adequate for simulating polymer solutions and colloidal suspensions, which are 
always over-damped, but may not be suitable for simulating complex plasmas, which are typically 



only lightly damped [|9j, 



it may not t 



LID. 



Currently, the most advanced BD algorithms [27, 



(to the best of our knowledge) are 
based on operator expansions of time pr opagato rs for Fokker-Plank or Kramers equations using 



Trotter decomposition technique IU5L 



16, 



L7L 



21. 



23 



27, 



28 



3 BOB. 



When 7 — > 0, these 



methods reduce to the class of so-called symplectic integrators for Hamiltonian dynamics, which 
have the unique properties of conserving the phase-space volume and being time reversible. These 
properties result in excellent stability properties and only very small energy drift on the long time 



23, 



27 



28 



29 



320 also use the standard Runge- 



Kutta algorithm to solve the deterministic part of dynamics and they, too, need multiple evaluations 
of deterministic forces per time step. As a consequence, they can become computationally quite 
expensive when used in simulation of many-particle systems, which restricts their applicability to 
systems of limited size, thus justifying further quest for more efficient algorithms of a comparable 
order. 
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On the other hand, the higher-order Predictor-Corrector (PC) methods, apparently, have not 
attracted much attention in BD simulation in comparison to their popularity in MD simulation 
J4, 14]. A PC method for BD was proposed by Ottinger in 1996 |5], and was subsequently adopted 
by Hutter [(]] in a simulation of colloidal suspensions. However, in their method the first step (pre- 
dictor) was just a first-order Euler-like scheme, which obviously does not utilize the full spectrum 
of advantages offered by the PC method. 

Consequently, we propose in the present paper a class of general purpose PC algorithms, which 
cover a wide range of friction coefficient 7 with higher-order accuracy, excellent stability, and 
very small energy drift on the long time scales. These methods reduce to Gear's methods for MD 



simulation in the limit of 7 — > [|4j, 



1411 . We shall discuss here only Gear- like algorithms that go 



up to the fifth-order, but extensions to higher orders should be quite straightforward, if required. 
The proposed Gear-like algorithms will be tested on simulation of a simple, one-particle system, 
and numerical comparisons will be made mainly with the Euler-like and Beeman-like (Verlet-like) 
methods because all these methods are based on a Taylor expansion of the deterministic force, 
Eq. (0, and they all require only one force evaluation per time step. As such, all these methods 
are specifically designed for simulations of large many-particle systems. Further discussion of ad- 
vantages of the proposed algorithms over those involving higher-order Runge-Kutta-like methods 



[21, 



23. 



27, 



28, 



3011 will be provided latter in the text. 



The paper is organized as follows. In Sec. II, we give detailed algorithms for conducting BD 
simulation. We next perform in Sec. Ill numerical tests for different methods (Euler-, Beeman- 



3j. 



and Gear-like) by using a stochastically damped harmonic oscillator (SDHO) model [33, 
Discussions of possible derivatives of the proposed BD method, as well as comparisons with the 
high-order Runge-Kutta-like methods are given in Sec. IV, which is followed by a brief conclusion 
in Sec. V. 



II. ALGORITHMS FOR BD SIMULATION 



A. General formula 



There are various ways to derive formulae for conducting BD simulation starting from 



Langevin equation Eq. dB. We shall follow here the strategy adopted in Refs. 



more recently in Refs. [[34 , 



and 



35|], because it is simple and straightforward, especially for readers 
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with some background in simulation but lacking background in SDEs. Following the technique of 



integrating the Langevin equation which is elaborated in Ref. 113411 . one realizes that Eq. (OQ) can be 
integrated exactly over a short time, based on certain assumptions about the deterministic accelera- 
tion a(t), thus giving updating formulae for BD simulation. The resulting formulae emphasize the 
fact that, under the assumptions that the stochastic acceleration in Eq. (OQ) is a Gaussian white noise 
and that 7 is constant, the two dynamic variables, v(t) and r(t), are actually normally distributed 
random variables themselves [or Gaussian random variables (GRVs)]. Consequently, according 



to the Normal linear transform theorem 11341 13511 . v(t) and r(t) are completely determined by 



their instantaneous means (vectors) and variances (scalars), and can be expressed, respectively, as 
follows 



v(t) = mean{v(t)} + Vvar{v(t)}N v (0,l) 



r(t) = mean{r(t)} + ^v<x{r(t)} N r (0, 1). (3) 

Here, N(0, 1) is a short-hand notation for a random vector [N = {N x ,N y ,N z } in three- 
dimensions (3D) or N = {N x , N y } in two-dimensions (2D)], whose components are mutually 
independent, standard (or unit) normal distributions (i.e., having the mean and variance 1). The 
subscripts v and r in Eq. ® emphasize that the two random normal vectors are associated with 
the velocity and position, respectively, and that they are different but correlated, as will be shown 
below. In the following, for the sake of simplicity, we shall occasionally denote the means by 
angular brackets, (v) and (r), and use standard deviations a v = ^/var{v} and a r = ^varjr}, 
instead of variances. 

We emphasize the significance of Eq. © because it provides the general updating formula for 
BD simulation. Namely, with Eq. ([3]), the process of solving the stochastic differential equation Eq. 
(OQ) and obtaining the two random variables v(t) and r(t) becomes simply a matter of determining 
the two sets of deterministic quantities: the means and the variances of v(t) and r(i), as well as 
the covariance between v(t) and r(t). Most importantly, the variances and covariances can be 
shown to be independent of the form of deterministic acceleration a(t). For variances, we find 



34. 



35D: 




o r = VvarM = U —Z- 1 - 2 — + — — — , (4) 
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where ks is the Boltzmann constant and T the temperature of the surrounding medium. In passing, 
we note that, when a Brownian particle is in equilibrium with the medium, T also defines its 
kinetic temperature. However, a Brownian particle may have a kinetic temperature different from 
T, which opens the possibilit y o f simulating non-equilibrium processes by using BD, such as 
particle stopping in a medium fl 1 IN . 

Further, the only non-zero covariances are co\{v x ,x\ = cov{v y , y} = cov{t> z , z] = cov{v , r}, 



and are given by (note, cov{i;, r} = (vr) — (v) (r) ) Q34 

cov{i;, r } = t^- (1 - 2e-* + e" 2 ^) . 
mrfi 



35] 



(5) 



Since the velocity and position are not independent, but rather jointly distributed normal variables, 
the updating formulae Eq. © can be further written as 113411 

r(t) = (r)+6 1 N 1 (0,l) + 6 2 N 2 (0,l), 

v(t) = (v) + ( 7,N 1 (0,l), (6) 

with b\ = cav{v,r}/a v and b 2 = y 'af — b\, where the unit normals Nj.(0, 1) and N 2 (0, 1) are 
now, by design, statistically (component-wise) independent of each other, and may be generated, 
e.g., by the Box-Muller method |36|] on a computer. 

It is worthwhile pointing out that the updating formulae Eq. © can be interpreted as a two-step 
algorithm: first, one calculates the means of the velocity and the position at time t and, second, 
one adds random increments of the velocity and position. Similar idea may also be seen in the 
Fokker-Planck based methods, e.g., in Ref. Il28ll . This notion is very important in the sense that 
the first step is no more than evaluation of deterministic Newton's equations with the addition 
of an average friction force, so that, in principle, any algorithm suitable for MD simulation (for 
example, Verlet, Beeman and even multi-step PC algorithms \a, 14]) may be used here. Moreover, 
with the second step being essentially independent of the first step, it can be always completed at 
the end of the time step as a correction to the previous step. More importantly, in the above two 
steps, the inaccuracy, or error, enters only in the first step through assumptions made regarding the 
deterministic acceleration a(t) Ill4ll . whereas in the second step, random increments of the velocity 
and position are added with the weights obtained from Eqs. (|4I5|) . which are exact and independent 
of the assumptions employed in the first step. 

Now, the only problem that remains is to determine the means of v(t) and r(i), which are 
closely dependent on the assumptions made about a(£), and which result in different simulation 
methods. We provide details for the Gear- like methods in the following. 



B. Gear-like Predictor-Corrector method 



Given the form of a(£), defined by the manner of how is the Taylor series in Eq. © truncated, 
(v) and (r) can be obtained in a number of ways by integrating Eq. ©. For simplicity, we follow 



here Lemons' methodology B41 I35I1. in which Eq. © is reduced to a set of deterministic (ordinary) 
differential equations by taking the statistical averages of both sides, giving 
d(r) 



dt 
d(v) 



(v) 

- 7 <v)+a(t), (7) 



dt 

where the fact that the Brownian acceleration A(t) is a Gaussian white noise with the mean 
was used. It is also important to note that the second equation in © may be further simplified by 
letting (v) = e~ 7< (v)', so that 

= a«. (8) 

dt 

Eq. © can be solved exactly with the help of Eq. ® and, subject to the initial conditions r , 

in) 

v , a , a , a , a and at t = 0, it yields 

(r) = r + civ i + c 2 a t 2 + c 3 a t 3 + c 4 a t 4 H h c„a[ ) n ~ 2) t' 1 H , 

(v) = c v + cia t + c 2 a t 2 + c 3 a t 3 + c 4 a t 4 + ■ ■ • + c^'V + • ■ ■ . (9) 

It is interesting to note that the coefficients in Eq. © obey very simple recursive relations. For 
n = 0, we have c = exp (— jt), while for n > 1, we find 

1 



[n - ljtyt 



(n- l)\c 



n-l 



j t (c n t n ) = c n _ x t n -\ (10) 

As expected, these coefficients reduce to those of a standard Taylor series when 7 — > 0. 

To construct the Gear-like method for BD simulation, one also needs higher order derivatives of 
the velocity (or position) I37L In principle, that may be attempted directly, by taking successively 
time derivatives on both sides of the second equation in Eq. ©, but that approach turns out to be 
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awkward. Instead, by using the equivalent equation Eq. ([8]), one finds 



-ft 



<*<v>' 



d_ 

dt 
<$_ 
dt 2 
f_ 
dt 3 



dt 
dt 

dt 

dt 



a(t) — a + a t + T^o^ 2 " 
a(t) = a + a t+ y&ot 2 

a(t) = a + a t H 

a(t) = a H 



s a »' 3 + 



(id 



Eq. (fTT]) . together with Eq. (Q, constitutes a Taylor-like series which can be used directly to 
construct the Gear-like PC method, as follows. 

In a direct analogy with simulations in deterministic systems, the Gear- like PC methods for BD 
simulation also include three stages, namely, predicting, force evaluating, and correcting 
3711 . The difference is that, one has to add random displacements of the velocity and position by 
using Eq. © at the end of a time step to complete the BD simulation. The basic procedure goes 
as follows. 

Predicting: In this stage, one has: 

r p = r + civ t + c 2 a t 2 + c 3 a t 3 + c 4 a t 4 + c 5 a t 5 
v p = c v + cia t + c 2 a t 2 + c 3 a t 3 + c 4 a i 4 



ag + a.gt + — a.gt 2 + — a'ot 3 



1 

3! 



a p 



a + SL t - 
a + a t 
a o? 



1 ... 2 

— r ant 
2! 



(12) 



where the superscript P is to indicate that these are quantities in the predicting stage. (For simplic- 
ity, we have dropped derivatives of a(t) higher than the third order, but extensions to higher orders 
are quite straightforward.) One may notice in Eq. (fT2l) that we have used in Eq. © the means of 
the position and velocity instead of using Taylor series of the position and velocity, as is normally 



done in the Gear method for MD simulation [|4 , 



1411 . Other than that, the remaining parts of this 



stage (derivatives of the force) are essentially the same as those in the MD simulation. 

Force evaluating: In this stage, the predicted position r p {t) is used to obtain a new force, that 
is, new acceleration a(t), and a difference between the predicted acceleration a p (t) and the new 
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acceleration a(t) is formed: 
Aa = [a(t) - a p (t)] 



(13) 



It can be seen that this step is exactly the same as in MD simulation 1141. 1 1411 . 

Correcting: In this stage, the above difference term is used to correct all predicted positions 
and their "derivatives", as follows: 



with 



v c t 
aFt 2 

2! 
a c £ 3 
~3T 
a c t 4 

4! 
a c t 5 
~5T 



AR2 



r p + 2c 2 a AR2 

v p t + c 1 a 1 AR2 
a, p t 2 

+ a 2 AR2 
+ a 3 AR2 
+ a 4 AR2 
+ a 5 AR2 



k p t 3 

3! 
z p t 4 

4! 

a p t 5 



5! 



Aat 2 



2! 



and the coefficients given in the following table: 





3th- order 


4th-order 


5th- order 


«0 


i 

6 


19 
120 


3 
16 




5 
6 


3 
4 


251 
360 


<y 2 


1 


1 


1 




1 
3 


1 

2 


11 
18 


a4 





1 

12 


1 

6 


«5 








1 

60 



This table is simply a reproduction of those in Refs. [|4j, [bj, |37|], and is provided here for com 



(14) 



(15) 



pleteness. By usingparameters in different columns, one may realize 3rd-, 4th-, and 5th-order (or 



4-, 5- and 6-value) |4j, LL4J Gear- like algorithms for BD simulation. Note that the first two lines in 
Eq. (fl~4~l) are slightly different from those in MD simulation, in order to restore the damping effect 
on deterministic acceleration, and to keep the consistency with the corresponding terms in Eq. © 
[or the first two lines in Eq. (fT2l) ]. 
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Adding random displacements: The above three stages are the same as those of Gear's algo- 
rithm for MD simulation. However, to complete the BD simulation, we have to use the updating 
formulae Eq. © to add random displacements of the velocity and position. It should be noted 
that now the corrected values v c and v c should be used in the places of (r) and (v) in Eq. ©, 
respectively. Also note that, to implement the simulation, one needs the initial conditions, r , v , 
a , a , ao and a , at t = 0. This poses a problem for the very first several steps of simulation, 
because a , a and a are undefined. The simplest way to get around this issue is to set all of them 
to be zero at the very first step, whereas their values will be updated during subsequent iterations. 
A better way would be to start the simulation by using a Runge-Kutta procedure for the first few 



steps 111411 . However, neither of these solutions of the problem will have any effect on the results 
in real, many-particle simulations. 

The above are the basic procedures for the Gear-like PC method for BD simulation. It should be 
restated that these formulae and the simulation stages are quite similar to those in the Gear method 



for MD simulation of Newton's equations [|4, 



14|], except for the expressions for the velocity and 



position in Eq. (fT2j) and Eq. (fl4j) . and for the addition of random displacement at the end of every 
time step. When 7 — > 0, the Gear-like method goes over to the Gear method for MD simulation. 
In the following we shall show that the performance of the proposed Gear-like method also shares 
many features with its counterpart in MD simulations. 

ffl. TESTING THE ALGORITHM 

In this section, we present some simple computational examples as testing cases for our Gear- 
like algorithm, and for the previously designed Euler-like and Beeman-like (Verlet-like) algo- 
rithms. For simplicity, we shall occasionally denote the Euler-like, Beeman-like, Verlet-like and 
Gear-like methods by EL, BL, VL, and GL, respectively. 



340, in which Eq. © may be 



We employ the model of a one-dimensional (ID) SDHO [|33L 
simplified to 

-T-r = v. (16) 
dt 



This model has been studied extensively and its solutions are well known Q33L 



34D. We have listed 



Chandrasekhar's results in the Appendix A to facilitate subsequent discussion. 
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FIG. 1 : (Color online) Deviations in the position, velocity and energy from their respective exact solutions 
within the first 100 time units. The left column is for 7 = 0.01 and the right one for 7 = 0.1. Here, GL-3, 
GL-4 and GL-5 denote the 3rd-, 4th-, and 5th-order Gear-like methods, respectively. The time step is fixed 
at 0.01. The (red) dashed lines, (blue) solid lines, and (black) solid lines show the deviations of energy, 
velocity and position, respectively. Note that the energy deviation is defined relative to the exact values by 
E(t)/E e {t) - 1. 

According to our discussion above, especially following Eq. ©, the task of a BD simulation is 
simply to predict the position r(t) and velocity v(t) of a Brownian particle at time t, given the set 
of initial conditions at time 0. As previously pointed out, r(t) and v(t) are normally distributed 
random vectors, whose variances and covariances are exact (at the level of Langevin equation), 
whereas their means are calculated numerically according to certain approximation schemes, such 
as EL, BL, VL or GL methods. From this point of view, the errors come only from evaluation of 
the means. In other words, the accuracy in calculating the means will decide the main performance 
of a BD simulation, so we begin our tests by considering the means. 
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FIG. 2: (Color online) Amplitudes of time-dependent deviations (or the largest deviations) in the position 
and velocity within 1000 time units for different methods. Here, EL and BL denote the Euler-like and 
Beeman-like methods, respectively. 7 = 0.01 and the time step of 0.05 are used in these tests. 



As usual |4|,ll4j], we firstly carry out simulations over certain periods and record, as a measure 
of accuracy, the deviations of r and v from their exact solutions without the inclusion of random 
displacements. We shall also judge the behavior of simulations by monitoring the quantity 



E(t) = v\t) + ujy(t) 



(17) 



which is proportional to the total energy of the Brownian particle. Note that, because of the exis- 
tence of damping, this energy is no longer a conserved quantity, but rather decreases exponentially 
with the factor exp(— 2^t). This may be seen from the exact solutions for the mean position and 
velocity listed in the Appendix C. In order to judge the energy conservation performance of numer- 
ical methods, we normalize E(t) by its exact counterpart E e (t), which is simply obtained from Eq. 
(fTTT) by substituting the exact values for the position and velocity. The resultant ratio E{t)/E e {t) 
should be a conserved quantity with the expected value of unity. 
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Without loss of generality, in all ensuing simulations we set k B T — 1, m — 1 and u = y/2/2, 
and choose the initial conditions to be r = 1 and vq = 0. These parameters are also used in the 
exact results based on Appendix A. 

The results of simulations for deviations in the position, velocity and energy from their exact 
solutions) are shown in Fig.Q]for the 3rd-, 4th-, and 5th-order Gear-like methods (denoted as GL- 
3, GL-4 and GL-5, respectively) within the first 100 time units. Note that the period of the damped 
oscillator is approximately 8.9 for both 7 = 0.01 and 7 = 0.1, and that the highest amplitude of 
velocity is approximately y/2/2. The oscillatory patterns of the deviations in position and velocity 
resemble approximately those of the exact solutions, but with much smaller amplitudes. The 
patterns of the energy deviations have a higher frequency (doubled) and seldom oscillate around 
0; rather, their oscillating centers drift in an approximately linear manner (in a semi-logarithmic 
plot). The slope of this line is usually defined as the energy drift, which is the main criterion testing 



the energy conservation in an algorithm |4|, [14J]. Of course, one would expect the slope to be as 
small as possible. Another criterion for the energy conservation is the amplitude of oscillations in 
the energy deviation, which is called the noise of the energy drift. Needless to say, the noise should 
be small, too, in accurate simulations. It is interesting to notice that some of the slopes in Fig. \T\ 
are negative, indicating that the energy is damped. This is quite different from the usual situation 
in the MD simulation where damping is absent by definition, and where the slopes of energy drift 
are always found positive. This, and other aspects of energy drift will be further discussed below. 

The amplitudes of the time-dependent deviations in position and velocity are shown in Fig. [2] 
for an extended time period of 1000 time units. One can easily see large differences in the ampli- 
tudes among different methods, with the patterns of the Euler-like and Beeman-like (Verlet-like) 
methods looking similar to those of the Gear-like method, but exhibiting orders of magnitude 
larger amplitudes than the Gear-like method. This is to be expected as, generally, all Gear-like 
methods have much better accuracy than Euler-like and Beeman-like (Verlet-like) methods. It 
should be furthermore pointed out that, except for deviations of the Euler-like method, which 
grow monotonously, all other deviations pass through broad maxima around the time of 200 time 
units. This indicates that the Euler-like method is not stable, while the others are. This illus- 
trates yet another aspect in which BD algorithms differ from those of MD simulation, where the 
deviations invariably grow with time |4J]. Namely, it is obvious that the finite damping 7 actu- 
ally stabilizes numerical algorithms (except for the Euler-like method). This is not surprising at 
all, because previous studies IU5L 1 171 . 13811 have demonstrated a possibility of stabilizing the MD 
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FIG. 3: (Color online) The largest deviation in the position from exact solution vs. the time step size in the 
first 20 time units, for different methods and different 7s. 

simulation by introducing a small, but finite damping. A similar idea is also well-known in the 



computational fluid dynamics where, for examp 



in fluid simulations to stabilize the computation [39] 



e, a small artificial viscosity is often introduced 



Next, we make some comparisons involving different time steps, as that is always the key 
issue in both the MD and BD simulations. Fig. [3] shows the largest deviations (defined by the 
largest deviation in the first 20 units) in position versus the size 5t of the time step for different 
methods and different friction coefficients 7. This figure is plotted in a double logarithmic scale, 
and the curves are nearly straight lines. The slope of these lines is called the apparent order [11411 - 
illustrating the dependence of error on the time step size. Namely, if the error is found to be 
proportional to 5t p , then the exponent p is the apparent order. It is found that for very small 7, for 
example 7 = 0.01, as shown in Fig. [3ja), the apparent orders are p E L ~ 1,Pbl ~ 2,Pgl-3 ~ 3.5, 
Pgl-4 ~ 4.2 and Pgl-5 ~ 4.6, respectively for the Euler-like, Beeman-like, 3rd-, 4th-, and 5th- 
order Gear-like methods. These values are very close to those from the MD simulations where 
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FIG. 4: (Color online) Time-dependent energy drift in GL-4 and GL-5 for different time step sizes St, with 
7 = 0.01 being fixed. 



damping is absent ||14|1 . When 7 increases, the absolute value of the error for Euler-like and 
Beeman-like methods decreases, while pel and Pel remain almost unchanged. On the other hand, 
Pgl-3, Pgl-4 and Pgl-5 slightly decrease with increasing 7. For example, for 7 = 1, as shown in 
Fig. He), Gear-like methods have the worst performance in accuracy: the apparent orders become 
now Pgl-3 ~ 3.1, Pgl-a ~ 3.9 and Pgl-5 ~ 4.0), but the amplitudes of their errors are still much 
smaller than those of the Beeman-like and Euler-like methods. 

We further provide a more detailed analysis of the energy conservation performances. We have 
already discussed the energy drift and noise for Gear-like methods in Fig.Q] In particular, negative 
energy drifts were noticed there for the 4th- and 5th-order Gear-like methods. However, this is 
not the case for all time step sizes, as illustrated in Fig. HI showing the energy drifts in the 4th- 
and 5th-order methods for different time step-sizes. We see that, when the step-size increases, 
the energy drift for GL-4 becomes increasingly positive, while that of GL-5 grows increasingly 
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FIG. 5: (Color online) Energy drift in one period of oscillation vs. the time step size for different methods 
and for 7 = 0.01 and 0.1. 

negative, until a very large time step is used (e.g., St = 0.6, not shown in the figure). On the other 
hand, the noise in the energy drift is always seen to increase with St. 

However, for GL-3, the trends in the energy drift are seen in Fig. [T] to be quite different, es- 
pecially for smaller 7. Namely, the oscillations in the energy deviation seem rather symmetrical 
about zero, with the amplitude of the noise being at least two orders of magnitude larger than the 
drift, as can be seen in the top left panel of Fig.[Q It is found that the Beeman-like method exhibits 
similar behavior as GL-3, which makes it very hard to do the least-squares fitting and linear re- 
gression analysis 111411 . On the other hand, for large 7, e.g. 7 > 1, the whole energy of the system 
is practically damped to zero in a few time units, so that it does not make much sense to discuss 
the energy drift in such cases. Therefore, in the following we shall omit the energy drift analysis 
for the BL and GL-3 methods altogether, as well as for other methods in the cases of large 7. 

Figured shows the energy drift over one period of oscillations versus the time step size | J 14f | for 
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FIG. 6: (Color online) Noise in the energy drift (defined by the average amplitude of oscillations in the first 
100 time units) vs. the time step size, for different methods and for 7 = 0.01 and 0.1. 

the EL, GL-4 and GL-5 methods with two different friction coefficients 7. Similarly, Fig. [6] shows 
the noise in the energy deviation (defined by the averaged amplitudes of oscillations in the first 100 
time units) versus the time step size, for these same methods with the same friction coefficients 
7. Since those figures are also plotted in a double logarithmic scale, we see that both the energy 
drift and the noise increase almost linearly with the time step size, and that the slopes of these 
curves are very close to the corresponding apparent orders, shown in Fig. [3] Moreover, the energy 
drifts for 7 = 0.1 in the GL-4 and GL-5 methods seem worse off than the corresponding drifts for 
7 = 0.01, but the magnitudes of the energy drift and noise are still very small. 

In the above, we have discussed the short-time performances of various BD methods in terms of 
accuracy and energy drift, without considering random displacements in the position and velocity. 
In the following, we perform several long-time tests including the random terms, in order to test 
the statistical properties of the simulation results. 
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FIG. 7: (Color online) Solid lines show the exact probability distribution functions (PDFs) for velocity and 
position given by Eq. (|A2|) . while symbols (circles and dots) show the statistics of 80, 000 samples taken 
from 800, 000 time units of simulation using the G-5 method. (Note that the PDFs have been scaled up, 
and their peaks are now unity.) Lower panels show the deviations of PDFs obtained in simulations from 
the exact ones, for the BL, GL-3, GL-4, and GL-5 methods. 7 = 0.01 and St = 0.1 are used in these 
simulations. 

In the long time limit, the oscillator approaches thermal equilibrium with its environment, so 
that both its velocity and position should become normally distributed. This may be easily seen 
from their exact solutions, Eq. (IA1I) . or more clearly from their long time limits, Eq. (IA2I) . which 
are plotted in Fig. [71 along with the statistical ensembles obtained from the simulation data. In this 
figure, solid lines show the probability distribution functions (PDFs) of the velocity and position 
given by Eq. (IA2I) . whereas symbols (circles and dots) show ensembles of 80, 000 sample points 
taken from 800, 000 time units of simulation. (Note that the first 2000 time units are used to 
allow the oscillator to reach an equilibrium, so that data from this period are not included in the 
ensembles shown in Fig.|7j). The dashed lines display deviations between the simulation data and 
the exact PDFs. We used 7 = 0.01 and St = 0.1 in simulations, allowing us to compare the results 
from the BL, GL-3, GL-4 and GL-5 methods, whereas the EL method is ruled out because it is not 
stable under these conditions. Since there are essentially no visible differences in the distribution 
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FIG. 8: (Color online) Deviations of PDFs for velocity and position for the EL method, with different time 
step sizes and damping coefficients. Note that the PDFs have been scaled up, and their peaks are now unity, 
as is shown in Fig. [7] 

functions generated by these four methods, we only show the PDF from the GL-5 method as an 
example, but display all deviations. We see that the deviations of all GL methods lie approximately 
within the range of (—3%, 3%), which is the statistical error for 80, 000 samples according to our 
experience. [We use the Box-Muller method to take 80, 000 samples, with the PDF given by Eq. 
(IA2I) . and compare the resultant PDF with the exact PDF. This gives a deviation in the range of 
(—3%, 3%).] The deviation from the BL method, as shown in the figure, is a bit larger: within 
the range of (—4%, 4%). With a smaller time step 5t and/or larger 7, all the above methods (EL, 
GL-3, GL-4, and GL-5) exhibit basically the same performance in the long time limit, and all their 
deviations are well within the range of the statistical error. 

Generally, a larger 7 would allow a larger time step 5t to be used for the same accuracy in 
the long time limit, and this is true for all methods including the EL method. For example, this 
method requires typically 7 > 1, along with the step size of about St = 0.1, enabling it to produce 
deviations in the PDF comparable to the other methods, as is shown in Fig. [8] 

IV. DISCUSSION 

A. Possible simplification 

We have provided detailed derivation of a class of Gear-like PC methods for performing BD 
simulation, and we tested them by using the SDHO model. In particular, we have compared 
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FIG. 9: (Color online) The correlation coefficient between the velocity and position vs. the time step size, 
for 7 = 0.01, 0.1, 1.0 and 10.0, respectively. 

the performances of different BD methods, namely, the Euler-like, Beeman-like (Verlet-like) and 
Gear-like methods, and discussed several similarities and differences between the BD and MD 
simulations. 

However, it may be sometimes computationally too demanding to implement a full BD simu- 
lation scheme, so we discuss briefly here some possible derivatives of the BD method. 

One possible simplification that might be made to the above BD simulation procedure is to 
neglect the correlation between the velocity and position, i.e., to neglect their covariance. This may 
be done by simply setting b\ = in Eq. ©, as many researchers actually do in their simulations. 
Therefore, we have also performed simulations without the covariance, and the results turned out 
essentially the same as those taking full account of the covariance. Physically, it is not easy to see 
why such a simplification should work at all. In fact, Fig. HI showing the correlation coefficient 
bi/a r versus the time step size for different 7, implies that the correlation between the velocity 
and position is not negligible at all. Therefore, it seems that monitoring the accuracy, energy 
conservation, and statistical properties of the ensembles of data from numerical results is not 
sufficient to fully evaluate performance of an algorithm. Further simulations with many-particle 
systems are currently being carried out to test if this simplification would affect any other physical 
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quantities of such systems. 

Yet another possible simplification to the above procedure is to apply a finite difference scheme 
for the deterministic force directly to Langevin equation (OQ), and numerically solve the resultant 
formulae like in MD simulation, but with the inclusion of corrections from the Langevin terms at 
the end of each time step. It can be proved that this simplification is equivalent to using linear 
approximations to the coefficients c n in the Taylor-like series of Eq. ©, for example, by letting 

Co w 1 - jt, 
1 

ci w 1 - -"ft, 

1 1 

Co ~ it, 

2 6 ' ' 

1 1 

(18) 

Clearly, when 7i <C 1, the above expressions would provide very good approximations to their 
real values given by Eq. ©. However, even if jt = 0.01, the use of this approximation would 
produce an error of magnitude around 10~ 4 , which is still much larger than the accuracy that the 
Gear-like methods can achieve. 



B. Comparison with High-order Runge-Kutta-like methods 

We provide here a brief discussion of the strengths and weaknesses of our algorithm in com- 
parison with some of the most recent developments in the area of high-order algorithms for BD 



28b, 



simulation, emphasizing that they typically involve Runge-Kutta methods H25L 

One should note that the algorithms based on decomposition of the exponential operator for 
time propagation provide superior stability of simulation, especially in the limit of low damping, 
owing to their unique properties of conserving the phase- space volume and being time reversible 
when simulating Hamiltonian dynamics problems. However, according to Bussi and Parrinello 



BHD, while such algorithms are derived with the aim of producing accurate trajectories up to a 
given order, they usually break down when a high damping is applied, which is the situation 
where time reversibility is irrelevant. Moreover, the design of such algorithms is not focused 
on the correctness of the ensemble generated in simulation 113 ID . which may become a serious 
shortcoming when simulating, e.g., equilibrium structure of a many-particle system interacting 
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with a bath. On the other hand, our method has demonstrated excellent performance in the long 
time limit, where equilibrium distribution functions were reproduced within the statistical error, 
both for the present single-particle model system, and for the many-particle applications in dusty 



plasmas 111, 



An additional advantage of the methods based on the time propagator with Fokker-Planck op- 
erator lies in the fact that various decomposition techniques make analysis quite transparent and 
allow for systematic strategies in developing higher-order schemes,the highest so far being the 
fourth-order numerical integrators developed by Drozdov and Brey Il27ll . and by Forbert and Chin 



H28H - However, it is precisely the decomposition of the exponentiated Fokker-Planck operator that 
seems to hinder further extensions to orders higher than fourth because of excessive growth in 
the complexity of such techniques with increasing order. On the other hand, simplicity of our 
algorithm makes extensions to orders higher than fifth quite straightforward, therefore offering an 
attractive pay-off for their lack of transparency. 

As mentioned in the Introduction, the most serious limitation of algorithms involving high- 



order Runge-Kutta-like methods, such as |25l 123, |28|], may arise when systems with large num 



bers of Brownian particles are simulated, where the need for force evaluations including all inter- 
particle interactions becomes the most critical issue in any simulation. For example, the method 
of Hershkovitz [250, when used in a ID simulation, requires four evaluations of force, four GRVs, 
and one evaluation of the force derivative per time step. Moreover, the method of Drozdov and 
Brey Q2VD , and the K4a method of Forbert and Chin Q28[| require three evaluations of force, one of 
the force derivative, and four GRVs for the former method while the latter method requires even 
eight GRVs. The K4b and K4c methods of Forbert and Chin 12811 require eight force evaluations 
and three GRVs per step. Yet the accuracy of these methods goes only up to the fourth order, 
whereas our method achieves easily the fifth-order accuracy with only one force evaluation and 
three GRVs per time step, which presents its greatest advantage for many-particle simulations. 

It should be mentioned, however, that the main weakness of our algorithm lies in the fact that 
it can handle only simple models of particle interactions with a bath, characterized with constant 
friction (diffusion) coefficients and with fluctuations driven by the Gaussian white noise. On the 



28D are 



other hand, the methods based on the time propagator decomposition techniques II27L 
capable of solving rather general SDEs, e.g., multivariable Langevin equation with configuration- 
dependent diffusion coefficients, or Kramers equation with colored noise. 
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V. CONCLUDING REMARKS 



In summary, a Gear-Like Predictor-Corrector algorithm was proposed for BD simulation. It 
has been tested by using a ID stochastically-damped-harmonic-oscillator model in terms of its 
accuracy, energy conservation, and the long-time statistical properties, and the results were com- 
pared with those obtained by using the Euler-like and Beeman-like (Verlet-like) methods. It was 
found that the present method exhibits much better performance in all the above tests. At the same 



time, when compared to the recent high-order Runge-Kutta-like methods 11251. 1271 128H. our algo- 
rithm promises superior efficiency for simulations in many-particle systems with constant friction 
(diffusion) coefficient and with Gaussian white noise, as only one force evaluation is needed per 
time step. 

We note that our Gear-like algorithm has already been used successfully for simulations of 
both equilibrium and nonequilibrium phenomena in strongly-coupled dusty plasmas, involving up 



mi pnei 

[lib 



to ~ 10 interacting dust particles II11L I12I1. More simulations with many-particle Yukawa systems 
are being carried out to further test the performance of our method, and the results will be reported 



elsewhere fl 1 3h 
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APPENDIX A: STOCHASTIC DAMPED HARMONIC OSCILLATOR 



The problem of a stochastic damped harmonic oscillator was first worked out in detail by Chan- 
drasekhar in 1943 113 311 . We provide here the list of his final results to facilitate the discussion in 
the main text. 
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where u)\ = a/7 2 — 4a; 2 ,. It should be noted that the above equations apply only when 7 > 2o; , 
i.e., in the over-damped case, while in the weakly-damped case, i.e., when 7 < 2u , one needs to 
simply redefine u>i as u)\ = a/4cJo — 7 2 , and replace cosh and sinh with cos and sin, respectively. 
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One can easily see from the above expressions that, when 7* 

r(t) = J^JV r (o,l) 

V muJ o 

v(t) = x n^N v (o,i), 

V m 



00, the solutions become: 



(A2) 
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